library(readxl)
library(car) # Para diagnósticos de regresión
library(MASS)10 Parcial 2 - Parte B
11 Introducción
Este documento detalla el proceso completo de ajuste y validación de un modelo de Regresión Lineal Múltiple. Se utiliza un conjunto de datos sobre los niveles de plasma en pájaros en función del tratamiento y el sexo para ilustrar la estimación de parámetros, pruebas de hipótesis, construcción de intervalos de confianza y un análisis de diagnóstico riguroso para identificar observaciones atípicas, de alta palanca e influyentes.
11.0.1 Librerías Requeridas
Se cargan los paquetes necesarios para el análisis.
11.0.2 Carga y Exploración de Datos
pajaros <- read_excel("pajaros.xlsx")
# View(pajaros)
plot(pajaros)
12 a) Estimación de los Parámetros del Modelo
Se ajusta un modelo de Regresión Lineal Múltiple de la forma: \[
\text{plasma}_i = \beta_0 + \beta_1 \cdot \text{treatment}_i + \beta_2 \cdot \text{sex}_i + \epsilon_i
\] Donde treatment y sex son variables binarias. Los parámetros (\(\beta_0, \beta_1, \beta_2\)) se estiman por Mínimos Cuadrados Ordinarios (MCO).
fit <- lm(plasma ~ 1 + treatment + sex, data = pajaros)
summary(fit)
Call:
lm(formula = plasma ~ 1 + treatment + sex, data = pajaros)
Residuals:
Min 1Q Median 3Q Max
-5.1037 -1.5087 -0.1075 1.0755 5.3564
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.5694 0.6156 23.667 < 2e-16 ***
treatment 17.0716 0.7108 24.016 < 2e-16 ***
sex -3.0342 0.7108 -4.269 0.000131 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.248 on 37 degrees of freedom
Multiple R-squared: 0.9415, Adjusted R-squared: 0.9383
F-statistic: 297.5 on 2 and 37 DF, p-value: < 2.2e-16# Intervalos de confianza para los coeficientes
confint(fit) 2.5 % 97.5 %
(Intercept) 13.322088 15.81676
treatment 15.631261 18.51186
sex -4.474525 -1.5939313 b) Prueba de Significancia Global del Modelo
Se evalúa si el modelo en su conjunto es significativo, es decir, si al menos uno de los predictores tiene una relación con la variable respuesta.
- \(H_0\): \(\beta_1 = \beta_2 = 0\) (El modelo no es significativo).
- \(H_1\): Al menos un \(\beta_j \neq 0\) (El modelo es significativo).
Esto se realiza mediante una prueba F.
# Opción 1: Usar el F-statistic y su p-valor de la salida de summary()
summary(fit)
Call:
lm(formula = plasma ~ 1 + treatment + sex, data = pajaros)
Residuals:
Min 1Q Median 3Q Max
-5.1037 -1.5087 -0.1075 1.0755 5.3564
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.5694 0.6156 23.667 < 2e-16 ***
treatment 17.0716 0.7108 24.016 < 2e-16 ***
sex -3.0342 0.7108 -4.269 0.000131 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.248 on 37 degrees of freedom
Multiple R-squared: 0.9415, Adjusted R-squared: 0.9383
F-statistic: 297.5 on 2 and 37 DF, p-value: < 2.2e-16# Opción 2: Comparar el modelo completo con un modelo reducido (solo intercepto) usando anova()
fit0 <- lm(plasma ~ 1, data = pajaros)
anova(fit0, fit, test = "F")Analysis of Variance Table
Model 1: plasma ~ 1
Model 2: plasma ~ 1 + treatment + sex
Res.Df RSS Df Sum of Sq F Pr(>F)
1 39 3193.4
2 37 187.0 2 3006.4 297.5 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusión: En ambas opciones, el p-valor es extremadamente pequeño (\(< 2.2e-16\)), por lo que se rechaza \(H_0\). El modelo es globalmente significativo.
14 d) Intervalo de Confianza para la Media
Se calcula un intervalo de confianza del 95% para el valor esperado de plasma (\(\mu_Y\)) para un pájaro con treatment=1 y sex=1.
La fórmula para el intervalo de confianza del valor esperado es: \[ \hat{\mu}_Y \pm t_{(1-\alpha/2, n-p)} \cdot EE(\hat{\mu}_Y) \] Donde \(EE(\hat{\mu}_Y) = \sqrt{\mathbf{x_0}^T (X^T X)^{-1} \mathbf{x_0} \hat{\sigma}^2}\).
x0_vector <- c(1, 1, 1) # Vector de condiciones: (intercepto, treatment=1, sex=1)
# Opción 1: Cálculo manual
n <- nrow(pajaros)
p <- length(coef(fit)) # Número de parámetros
mu_hat <- t(x0_vector) %*% as.matrix(fit$coeff)
ee <- sqrt(t(x0_vector) %*% vcov(fit) %*% x0_vector)
perct <- qt(0.975, df = n - p)
LI <- mu_hat - perct * ee
LS <- mu_hat + perct * ee
c(LI = LI, LS = LS) LI LS
27.35942 29.85409 # Opción 2: Usando la función predict()
nuevos_datos <- data.frame(treatment = 1, sex = 1)
predict(fit, nuevos_datos, interval = "confidence", level = 0.95) fit lwr upr
1 28.60675 27.35942 29.8540915 e) Prueba de Interacción entre treatment y sex
Se evalúa si el efecto del tratamiento sobre el plasma es diferente para machos y hembras, ajustando un modelo con un término de interacción.
- \(H_0\): No hay interacción (\(\beta_{treatment:sex} = 0\)).
- \(H_1\): Hay interacción (\(\beta_{treatment:sex} \neq 0\)).
# Se ajusta el modelo con interacción
fit_interaccion <- lm(plasma ~ 1 + treatment + sex + treatment:sex, data = pajaros)
summary(fit_interaccion)
Call:
lm(formula = plasma ~ 1 + treatment + sex + treatment:sex, data = pajaros)
Residuals:
Min 1Q Median 3Q Max
-4.6803 -1.1562 -0.3421 0.9375 5.7798
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.1460 0.7067 20.017 <2e-16 ***
treatment 17.9183 0.9994 17.929 <2e-16 ***
sex -2.1874 0.9994 -2.189 0.0352 *
treatment:sex -1.6936 1.4134 -1.198 0.2387
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.235 on 36 degrees of freedom
Multiple R-squared: 0.9437, Adjusted R-squared: 0.939
F-statistic: 201.1 on 3 and 36 DF, p-value: < 2.2e-16# Se comparan los modelos con y sin interacción usando una prueba F parcial
anova(fit, fit_interaccion, test = "F")Analysis of Variance Table
Model 1: plasma ~ 1 + treatment + sex
Model 2: plasma ~ 1 + treatment + sex + treatment:sex
Res.Df RSS Df Sum of Sq F Pr(>F)
1 37 186.96
2 36 179.79 1 7.1703 1.4358 0.2387Conclusión: El p-valor (0.898) es grande, por lo que no se rechaza \(H_0\). No hay evidencia de un efecto de interacción significativo.
16 f) Prueba de Hipótesis Lineal: \(\beta_1 + \beta_2 = 0\)
Se prueba si la suma de los coeficientes de treatment y sex es igual a cero. - \(H_0\): \(\beta_1 + \beta_2 = 0\). - \(H_1\): \(\beta_1 + \beta_2 \neq 0\).
Se utiliza un estadístico t basado en la combinación lineal de los coeficientes: \[ t_c = \frac{\mathbf{a}^T \hat{\boldsymbol{\beta}} - c}{\sqrt{\hat{\sigma}^2 \mathbf{a}^T (X^T X)^{-1} \mathbf{a}}} \] donde \(\mathbf{a}^T = [0, 1, 1]\) y \(c=0\).
a <- c(0, 1, 1) # Define la combinación lineal
n <- nrow(pajaros)
p <- length(coef(fit))
t_C <- (t(a) %*% as.matrix(fit$coeff)) / sqrt(t(a) %*% vcov(fit) %*% a)
# Cálculo del p-valor
pval_n <- 2 * (1 - pt(abs(t_C), df = n - p))
c(t_calculado = t_C, p_valor = pval_n) t_calculado p_valor
1.396362e+01 2.220446e-16 Conclusión: El p-valor es muy pequeño, por lo que se rechaza \(H_0\).
17 h) Análisis de Diagnóstico: Puntos Influyentes
Se realiza un análisis de diagnóstico para identificar observaciones que puedan tener un efecto desproporcionado en el ajuste del modelo.
17.1 Alta Palanca (Leverage)
Las observaciones con alto apalancamiento son aquellas con valores atípicos en los predictores (\(X\)). Un punto tiene alta palanca si su valor de hat, \(h_{ii}\), es mayor que \(2p/n\) o \(3p/n\).
n <- nrow(pajaros)
p <- length(coef(fit)) # Número de parámetros (b0, b1, b2)
# Gráfico de los valores de apalancamiento (hat values)
plot(hatvalues(fit), type = "h", main = "Índice de Valores de Apalancamiento")
abline(h = 2 * p / n, col = "blue", lty = 2)
abline(h = 3 * p / n, col = "red", lty = 2)
legend("topright", legend = c("Corte 2p/n", "Corte 3p/n"), col = c("blue", "red"), lty = 2)
# Identificar los puntos con mayor apalancamiento
head(sort(hatvalues(fit), decreasing = TRUE)) 3 4 5 6 7 8
0.075 0.075 0.075 0.075 0.075 0.075 17.2 Observaciones Atípicas (Outliers)
Los outliers son observaciones con un gran error de predicción (residual grande). Se identifican usando los residuales estudentizados. Una observación se considera atípica si \(|r_{i,stud}| > 3\).
# Residuales estudentizados
stud_res <- studres(fit)
head(sort(abs(stud_res), decreasing = TRUE)) 2 21 34 14 3 19
2.675926 2.553484 2.526613 1.792766 1.378081 1.340180 # Visualización con Boxplot
boxplot(stud_res, main = "Boxplot de Residuales Estudentizados")
Conclusión: No se observan residuales estudentizados con un valor absoluto mayor a 3, por lo que no hay outliers claros.
17.3 Observaciones Influyentes
Una observación influyente es aquella que, si se elimina, cambia significativamente el ajuste del modelo. Combina tanto el apalancamiento como el tamaño del residual.
- Distancia de Cook: Mide el efecto de eliminar la i-ésima observación. Una regla general es que \(D_i > 4/(n-p)\) es potencialmente influyente.
- Gráfico de Influencia: Visualiza simultáneamente el apalancamiento, el residual estudentizado y la Distancia de Cook.
par(mfrow = c(1, 2))
# Gráfico de la Distancia de Cook
corte <- 4 / (n - p)
plot(fit, which = 4, cook.levels = corte)
abline(h = corte, lty = 2, col = "red")
# Gráfico de Influencia (burbujas)
influencePlot(fit, id.method = "identify", main = "Gráfico de Influencia", sub = "El tamaño del círculo es proporcional a la D_Cook")Warning in plot.window(...): "id.method" is not a graphical parameterWarning in plot.xy(xy, type, ...): "id.method" is not a graphical parameterWarning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
a graphical parameter
Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
a graphical parameterWarning in box(...): "id.method" is not a graphical parameterWarning in title(...): "id.method" is not a graphical parameterWarning in plot.xy(xy.coords(x, y), type = type, ...): "id.method" is not a
graphical parameter
StudRes Hat CookD
2 2.6759260 0.075 0.165905553
3 -1.3780809 0.075 0.050109542
4 -0.3778034 0.075 0.003949214
21 2.5534838 0.075 0.15334516818 Comparación de modelos con y sin puntos influyentes
Se comparan los coeficientes del modelo original con un modelo ajustado sin las observaciones identificadas como más influyentes.
# Puntos identificados como influyentes (ejemplo)
puntos_influyentes <- c(2, 21, 34)
# Modelo sin las observaciones influyentes
fit_sin_influyentes <- update(fit, subset = -puntos_influyentes)
# Comparación de resúmenes
summary(fit)
Call:
lm(formula = plasma ~ 1 + treatment + sex, data = pajaros)
Residuals:
Min 1Q Median 3Q Max
-5.1037 -1.5087 -0.1075 1.0755 5.3564
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.5694 0.6156 23.667 < 2e-16 ***
treatment 17.0716 0.7108 24.016 < 2e-16 ***
sex -3.0342 0.7108 -4.269 0.000131 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.248 on 37 degrees of freedom
Multiple R-squared: 0.9415, Adjusted R-squared: 0.9383
F-statistic: 297.5 on 2 and 37 DF, p-value: < 2.2e-16summary(fit_sin_influyentes)
Call:
lm(formula = plasma ~ 1 + treatment + sex, data = pajaros, subset = -puntos_influyentes)
Residuals:
Min 1Q Median 3Q Max
-2.9710 -0.9354 -0.0917 1.1142 3.6393
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.8348 0.4864 28.444 < 2e-16 ***
treatment 17.3736 0.5568 31.204 < 2e-16 ***
sex -2.1740 0.5568 -3.905 0.000425 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.692 on 34 degrees of freedom
Multiple R-squared: 0.967, Adjusted R-squared: 0.965
F-statistic: 498 on 2 and 34 DF, p-value: < 2.2e-16